Upregulated Expression of ERBB2/HER2 in Multiple Myeloma as a Predictor of Poor Survival Outcomes

The main goal of the present study was to examine if the RNA-sequencing (RNAseq)-based ERBB2/HER2 expression level in malignant plasma cells from multiple myeloma (MM) patients has clinical significance for treatment outcomes and survival. We examined the relationship between the RNAseq-based ERBB2 messenger ribonucleic acid (mRNA) levels in malignant plasma cells and survival outcomes in 787 MM patients treated on contemporary standard regimens. ERBB2 was expressed at significantly higher levels than ERBB1 as well as ERBB3 across all three stages of the disease. Upregulated expression of ERBB2 mRNA in MM cells was correlated with amplified expression of mRNAs for transcription factors (TF) that recognize the ERBB2 gene promoter sites. Patients with higher levels of ERBB2 mRNA in their malignant plasma cells experienced significantly increased cancer mortality, shorter progression-free survival, and worse overall survival than other patients. The adverse impact of high ERBB2 expression on patient survival outcomes remained significant in multivariate Cox proportional hazards models that accounted for the effects of other prognostic factors. To the best of our knowledge, this is the first demonstration of an adverse prognostic impact of high-level ERBB2 expression in MM patients. Our results encourage further evaluation of the prognostic significance of high-level ERBB2 mRNA expression and the clinical potential of ERBB2-targeting therapeutics as personalized medicines to overcome cancer drug resistance in high-risk as well as relapsed/refractory MM.

While the expression and function of the epidermal growth factor receptor (EGFR)/ERBB1 has been extensively studied in normal and malignant cells of epithelial and mesenchymal origins, very little is understood regarding its expression in normal or malignant lymphoid cells. However, recent studies indicated that lymphoid cells also express ERBB1 [16][17][18][19][20]. Notably, ERBB1 ligands AREG and EGF have been shown to stimulate the proliferation of multiple myeloma (MM) cells [21,22], and inhibiting ERBB1 caused cytotoxicity to MM cells [23,24]. Our recent studies demonstrated that ERBB1 expression is upregulated in MM cells, and higher levels of ERBB1 are associated with poor treatment outcomes and survival in newly diagnosed MM [25].
Although several ligands have been identified for ERBB1, ERBB3, and ERBB4, a specific ligand for the receptor kinase ERBB2 has not been identified [2,4]. However, ERBB2/HER2 serves as a co-receptor for other ERBB proteins that is capable of forming heterodimers and thereby facilitating the generation of an amplified intracellular biochemical signal following binding of a specific ligand to the extracellular domain of its heterodimerization partner [4]. For example, ERBB1-binding ligands, such as EGF, amphiregulin (AR), and transforming growth factor-α (TGFα) trigger a potent signal after binding to the ERBB1 portion of ERBB1×ERBB2 heterodimers [4]. The primary purpose of the present bioinformatics study was to examine if the RNA-sequencingbased ERBB2/HER2 expression level in malignant plasma cells from MM patients has clinical significance for treatment outcomes and survival. To this end, we used the Multiple Myeloma Research Foundation (MMRF)-CoMMpass RNA sequencing (RNAseq) dataset generated in patients treated on contemporary standard regimens. Upregulated expression of ERBB2 in MM cells was correlated with amplified expression of transcription factors (TF) that recognize the ERBB2 gene promoter sites. Notably, high-level ERBB2/HER2 expression was associated with increased cancer mortality and significantly worse overall survival (OS) of MM patients.

ERBB2/HER2 mRNA Is Expressed at Significantly Higher Levels Than EGRF/ERBB1 mRNA and ERBB3 mRNA across All Clinical Stages of MM
The ERBB2/HER2 gene was abundantly expressed in malignant plasma cells and none of the 766 MM patients with ISS staging information had zero alignments to the ERBB2/HER2 gene ( Figure 1A). No significant differences in ERBB2/HER2 expression levels were observed in the pairwise comparisons between the three stages ( Figure 1B). Notably, the levels for ERBB2/HER2 mRNA in the pooled set of 766 MM patients (Mean ± SEM = 7.71 ± 0.03) was significantly higher than the levels for ERBB1 mRNA (Mean ± SEM = 4.08 ± 0.04; t-test, p < 0.0001) or ERBB3 mRNA (Mean ± SEM = 4.09 ± 0.03; p < 0.0001).

Amplified ERBB2/HER2 Expression in Malignant Plasma Cells from MM Patients Is Associated with Poor PFS Outcomes
We first evaluated the potential impact of high-level ERBB2/HER2 mRNA expression on PFS outcome in 669 evaluable MM patients ( Figure S1A). ERBB2 high patients (i.e., 335 patients with the top 50% of the highest observed expression level for FPKM-UQ aligned to the ERBB2 sequence) had significantly shorter median PFS than ERBB2 low patients (i.e., 334 patients with the bottom 50% of the expression level for FPKM-UQ aligned to the ERBB2 sequence) (27 months vs. 36 months, log-rank chi-square = 4.85, p-value = 0.028). We confirmed the adverse impact of high ERBB2/HER2 expression on PFS using both univariate and multivariate Cox regression models for 647 evaluable patients ( Figure 3A). Further, in the multivariate Cox regression model, ERBB2 mRNA expression level as a linear covariate exhibited a significant 1.54-fold increase in the hazard ratio for each unit increase in FPKM-UQ counts (p < 0.001) that persisted with the inclusion of other prognostic factors, including ISS stage, age, as well as serum albumin and beta 2 microglobulin levels ( Figure S1B).

Amplified ERBB2/HER2 Expression in Malignant Plasma Cells from MM Patients Is Associated with Poor PFS Outcomes
We first evaluated the potential impact of high-level ERBB2/HER2 mRNA expression on PFS outcome in 669 evaluable MM patients ( Figure S1A). ERBB2 high patients (i.e., 335 patients with the top 50% of the highest observed expression level for FPKM-UQ aligned to the ERBB2 sequence) had significantly shorter median PFS than ERBB2 low patients (i.e., 334 patients with the bottom 50% of the expression level for FPKM-UQ aligned to the ERBB2 sequence) (27 months vs. 36 months, log-rank chi-square = 4.85, p-value = 0.028). We confirmed the adverse impact of high ERBB2/HER2 expression on PFS using both univariate and multivariate Cox regression models for 647 evaluable patients ( Figure 3A). Further, in the multivariate Cox regression model, ERBB2 mRNA expression level as a linear covariate exhibited a significant 1.54-fold increase in the hazard ratio for each unit increase in FPKM-UQ counts (p < 0.001) that persisted with the inclusion of other prognostic factors, including ISS stage, age, as well as serum albumin and beta 2 microglobulin levels ( Figure S1B).

Amplified ERBB2/HER2 Expression in Malignant Plasma Cells from MM Patients Is Associated with Poor OS
We next sought to determine the effect of amplified ERBB2/HER2 mRNA expression on OS in 787 MM patients considering all deaths from any cause in our analysis. ERBB2 high patients (i.e., 394 patients with the top 50% of the highest observed expression level for FPKM-UQ aligned to the ERBB2 sequence) had a significantly shorter OS than ERBB2 low patients (i.e., 393 patients with the bottom 50% of the observed expression level for FPKM-UQ aligned to the ERBB2 sequence) had significantly worse OS outcomes. The cumulative proportion of ERBB2 high patients remaining alive at 60 months was 50.1 ± 6.7%, whereas the cumulative proportion of ERBB2 low patients remaining alive at 60 months was 69.5 ± 7.2% (log-rank chi-square = 8.82, p-value = 0.003) ( Figure 4A). We confirmed the adverse impact of high ERBB2 expression on OS using both univariate and multivariate Cox regression models for 767 evaluable patients ( Figure 3B). Further, in the multivariate Cox regression model, ERBB2 mRNA expression level as a linear covariate exhibited a significant two-fold increase in the hazard ratio for each unit increase in FPKM-UQ counts (p < 0.001) that persisted with the inclusion of other prognostic factors, including ISS stage, age, as well as serum albumin and beta 2 microglobulin levels ( Figure 4B).
We next assessed the OS of four distinct patient groups defined by co-expression levels of ERBB2 and ERBB1: (i) ERBB1 high /ERBB2 high (top 50th percentile of the highest observed expression levels for both ERBB1 and ERBB2 mRNA expression levels) (n = 195), (ii) ERBB1 low /ERBB2 low (bottom 50th percentile of the lowest observed expression levels for both ERBB1 and ERBB2 (n = 194), (iii) ERBB1 low /ERBB2 high (n = 199), and (iv) ERBB1 high /ERBB2 low (n = 199). ERBB1 low /ERBB2 low patients representing 24.7% of the 787 evaluable patients had significantly fewer deaths and consequently the best 5-year survival of 83.3% when compared to the remaining three groups of patients ( Figure S2).

Amplified ERBB2/HER2 Expression in Malignant Plasma Cells from MM Patients Is Associated with Increased Cancer-Related Mortality
We next evaluated the potential impact of high-level ERBB2/HER2 mRNA expression on cancer-related mortality (CRM) in 716 evaluable MM patients ( Figure S3A). ERBB2 high patients (i.e., 358 patients with the top 50% of the observed expression level for FPKM-UQ aligned to the ERBB2 sequence) had a significantly higher CRM than ERBB2 low patients (i.e., 358 patients with the bottom 50% of the observed expression level for FPKM-UQ aligned to the ERBB2 sequence) (37.3 ± 7.2% vs. 22.7 ± 8.5%; log-rank chi-square = 6.17, p-value = 0.013). We confirmed the adverse impact of high ERBB2 expression on CRM using a multivariate Cox regression model for 695 evaluable patients. ERBB2 mRNA expression level as a linear covariate exhibited a significant 2.23-fold increase in the hazard ratio for each unit increase in FPKM-UQ counts (p < 0.001) that persisted with the inclusion of other prognostic factors, including ISS stage, age, as well as serum albumin and beta 2 microglobulin levels ( Figure S3B). We next assessed the OS of four distinct patient groups defined by co-expression levels of ERBB2 and ERBB1: (i) ERBB1 high /ERBB2 high (top 50th percentile of the highest observed expression levels for both ERBB1 and ERBB2 mRNA expression levels) (n = 195), (ii) ERBB1 low /ERBB2 low (bottom 50th percentile of the lowest observed expression levels for both ERBB1 and ERBB2 (n = 194), (iii) ERBB1 low /ERBB2 high (n = 199), and (iv) ERBB1 high /ERBB2 low (n = 199). ERBB1 low /ERBB2 low patients representing 24.7% of the 787 evaluable patients had significantly fewer deaths and consequently the best 5-year survival of 83.3% when compared to the remaining three groups of patients ( Figure S2). Cumulative proportions surviving at 60 mo: 50.1 ± 6.7% for ERBB2 high and 69.5 ± 7.2% for ERBB2 low . (B) Depicted is the forest plot comparing the HRs for each covariate included in a multivariate Cox proportional hazards model, namely, ERBB2 mRNA level, ISS stage, age, as well as serum albumin and beta 2 microglobulin levels in 767 evaluable patients. ERBB2 mRNA expression level exhibited a significant 1.96-fold increase in the hazard ratio for each unit increase in FPKM-UQ counts (p < 0.001). Significant p-values are indicated by *, **, and *** for p < 0.05, p < 0.01, and p < 0.001, respectively. We next assessed the CRM of four distinct patient groups defined by co-expression levels of ERBB2 and ERBB1: (i) ERBB1 high /ERBB2 high (the top 50% of the observed expression levels for both ERBB1 and ERBB2 mRNA expression levels) (n = 181), (ii) ERBB1 low /ERBB2 low (the bottom 50% of the observed expression levels for both ERBB1 and ERBB2 (n = 181), (iii) ERBB1 low /ERBB2 high (n = 177), and (iv) ERBB1 high /ERBB2 low (n = 177). ERBB1 low /ERBB2 low patients representing 25.3% of the 716 evaluable patients had significantly fewer cancerrelated deaths and, consequently, the lowest 5-year CRM (7.2 ± 2.1%) when compared with the remaining three groups of patients (5-year CRM: ERBB1 high /ERBB2 high = 39.4 ± 10%, ERBB1 high /ERBB2 low = 38.1 ± 14.4%, and ERBB1 low /ERBB2 high = 33.8 ± 9%) ( Figure S4).

A STRING Model of HER2/ERBB2-Regulated Signaling Network Reveals 8 ERBB2-Associated Signaling Molecules as Poor OS Indicators
We compared the OS outcomes curves for 394 MM patients with high levels of mRNA expression for ERBB2 and associated genes coding for proteins associated with ERBB2 (i.e., patients with the top 50% of the observed expression level for FPKM-UQ values for each gene) with the OS outcomes of 393 MM patients with low mRNA expression levels for each gene (i.e., patients with the bottom 50% of the observed expression level for FPKM-UQ values). Depicted in Figure 5
ERBB2/HER2 serves as a co-receptor for ERBB1 and as a heterodimer with ERBB1 contributing to the generation of an amplified intracellular biochemical signal after binding of an ERBB1 ligand to the ERBB1 portion of the ERBB1×ERBB2 heterodimer [4]. Upregulated ERBB2/HER2 expression could contribute to the formation of ERBB1×ERBB2 heterodimers and thereby potentiate the ability of ERBB1 ligands in the MM TME, such as EGF and AR to promote the survival, proliferation, and dissemination of MM cells. Patients with higher levels of ERBB2 mRNA in their malignant plasma cells experienced significantly increased CRM, shorter PFS, and worse OS than other patients. The adverse impact of high ERBB2/HER2 expression on patient survival outcomes remained significant in multivariate Cox proportional hazards models that accounted for the effects of other prognostic factors.
In MM cells with low-level ERBB1 expression unlikely to form ERBB1xERBB1 homodimers abundantly, abundant expression of ERBB2 could compensate for the low-level ERBB1 expression by enabling the formation of ERBB1×ERBB2 heterodimers. It is note-worthy that the CRM of ERBB1 low ERBB2 high patients was as high as the CRM of ERBB1 high patients with low or high ERBB2 levels ( Figure S4). We propose a model according to which ERBB1 signals transmitted by activated ERBB1×ERBB1 homodimers or ERBB1×ERBB2 heterodimers contribute to the demonstrated adverse prognostic effects of ERBB1 [25] and ERBB2/HER2 (present study) in MM ( Figure 6). ERBB2 exhibited direct links with eight signaling molecules that were associated with significant increases in hazard ratios for poor OS (Table 1, Figure 5). In agreement with such a model, ERBB1 lowand ERBB2 low patients whose MM cells would be less likely to form biologically meaningful quantities or ERBB1xERBB1 homodimers or ERBB1×ERBB2 heterodimers had the lowest 5-year cancer-related mortality (7.2 ± 2.1%) when compared with the remaining patients ( Figure S4). inhibit both ERBB1 and ERBB2 and may have clinical potential for high-risk or R/R MM patients. Figure 6. Amplified HER2/ERBB2 expression and directly associated signaling molecules identified via the STRING algorithm resulting in poor overall survival in MM patients. The HER2/ERBB2 receptor (blue symbol) does not bind to the ligands but dimerizes with the EGFR/ERBB1 receptor (red symbol) to activate downstream signal transduction pathways. Dimers of EGFR/ERBB1 bind to EGFR/ERBB1 ligands (green circle; a select group of ligands shown in the green box including EGF, amphiregulin (AREG), transforming growth factor-alpha (TGFα), and heparin-binding EGF-like growth factor (HB-EGF)) to send signals into the cell. Depicted are the 8 signaling molecules directly associated with ERBB2 from the proposed interaction network identified via the STRING algorithm that also exhibited poor OS outcomes at high levels of mRNA expression.
In this study, we demonstrated that ERBB2/HER2 co-receptor mRNA was expressed in malignant plasma cells from MM patients at significantly higher levels than ERBB1 as well as ERBB3 across all three stages of the disease. We tested the hypothesis that the amplified ERBB2 mRNA levels in MM cells may be driven by increased expression of specific TFs that transcriptionally activate ERBB2 expression. Our findings demonstrated for the first time that increased ERBB2 mRNA expression shows statistically significant correlations with upregulated mRNA levels of several such TFs, including ETV, SP1, CE-BPB, PHF8, and TBP. We propose a model according to which the observed upregulation of ERB2/HER2 mRNA expression in MM cells is driven transcriptionally by several TFs ( Figure 2C). The presented data expand our current knowledge regarding the networks of signaling pathways affecting the biology and clinical outcome of MM as well as emerging new molecular targets for chemotherapy-drug-resistant MM [60][61][62][63][64][65][66][67][68][69].
Our hypothesis-generating study suffers from a number of limitations including its primary focus on bioinformatics-based analyses without integrated laboratory testing of ERBB2 mRNA levels using other methods in addition to RNAseq and lack of data to show that MM cells from ERBB2 high patients with high levels of ERBB2 mRNA also express ERBB2 protein at higher levels than ERBB2 low patients. Our results need to be verified in a prospective hypothesis-testing study with integrated RNAseq and biochemical testing in a large MM patient population. A multivariate analysis of the prognostic significance of ERBB2 overexpression in relationship to other prognostic factors and especially patient treatment will be very important.

Processing and Analysis of the Multiple Myeloma Research Foundation (MMRF)-CoMMpass RNA Sequencing (RNAseq) Dataset
The processing of the MMRF-CoMMpass RNAseq dataset was performed, as previously detailed [25]. Pre-processed RNAseq data files for 787 cases used in our analysis are Figure 6. Amplified HER2/ERBB2 expression and directly associated signaling molecules identified via the STRING algorithm resulting in poor overall survival in MM patients. The HER2/ERBB2 receptor (blue symbol) does not bind to the ligands but dimerizes with the EGFR/ERBB1 receptor (red symbol) to activate downstream signal transduction pathways. Dimers of EGFR/ERBB1 bind to EGFR/ERBB1 ligands (green circle; a select group of ligands shown in the green box including EGF, amphiregulin (AREG), transforming growth factor-alpha (TGFα), and heparin-binding EGF-like growth factor (HB-EGF)) to send signals into the cell. Depicted are the 8 signaling molecules directly associated with ERBB2 from the proposed interaction network identified via the STRING algorithm that also exhibited poor OS outcomes at high levels of mRNA expression. Several ERBB2/HER-2 targeting therapeutics have been approved for the treatment of solid tumors, such as breast cancer and lung cancer [57][58][59]. Targeting ERBB2 with FDA-approved small molecule ERBB2 inhibitors as well as monoclonal antibodies could potentially improve the treatment options for high-risk or relapsed/refractory (R/R) MM patients. Our results presented herein warrant further examination of the therapeutic potential of already FDA-approved small molecule inhibitors as well as monoclonal antibodies targeting ERBB2 that could be repurposed for use in high-risk or R/R MM. The evaluation of the clinical potential of ERBB2/HER-2 targeting therapeutics, including monoclonal antibodies trastuzumab and pertuzumab, antibody-drug conjugates trastuzumab emtansine and trastuzumab deruxtecan, and small molecule TKI tucatinib alone and in combination with standard anti-MM drugs, as well as ERBB1-targeting antibodies and TKI, would seem warranted. The dual function of TKI neratinib and palatinib inhibit both ERBB1 and ERBB2 and may have clinical potential for high-risk or R/R MM patients.
In this study, we demonstrated that ERBB2/HER2 co-receptor mRNA was expressed in malignant plasma cells from MM patients at significantly higher levels than ERBB1 as well as ERBB3 across all three stages of the disease. We tested the hypothesis that the amplified ERBB2 mRNA levels in MM cells may be driven by increased expression of specific TFs that transcriptionally activate ERBB2 expression. Our findings demonstrated for the first time that increased ERBB2 mRNA expression shows statistically significant correlations with upregulated mRNA levels of several such TFs, including ETV, SP1, CEBPB, PHF8, and TBP. We propose a model according to which the observed upregulation of ERB2/HER2 mRNA expression in MM cells is driven transcriptionally by several TFs (Figure 2C). The presented data expand our current knowledge regarding the networks of signaling pathways affecting the biology and clinical outcome of MM as well as emerging new molecular targets for chemotherapy-drug-resistant MM [60][61][62][63][64][65][66][67][68][69].
Our hypothesis-generating study suffers from a number of limitations including its primary focus on bioinformatics-based analyses without integrated laboratory testing of ERBB2 mRNA levels using other methods in addition to RNAseq and lack of data to show that MM cells from ERBB2 high patients with high levels of ERBB2 mRNA also express ERBB2 protein at higher levels than ERBB2 low patients. Our results need to be verified in a prospective hypothesis-testing study with integrated RNAseq and biochemical testing in a large MM patient population. A multivariate analysis of the prognostic significance of ERBB2 overexpression in relationship to other prognostic factors and especially patient treatment will be very important.

Processing and Analysis of the Multiple Myeloma Research Foundation (MMRF)-CoMMpass RNA Sequencing (RNAseq) Dataset
The processing of the MMRF-CoMMpass RNAseq dataset was performed, as previously detailed [25]. Pre-processed RNAseq data files for 787 cases used in our analysis are accessible via the GDC portal (https://gdc.cancer.gov/about-gdc/contributed-genomic-datacancer-research/foundation-medicine/multiple-myeloma-research-foundation-mmrf, accessed on 20 March 2023). The analytical path for RNAseq datasets employed standardized mRNA quantification techniques to enable meta-analysis across multiple projects, as previously described in detail [25]. We utilized the latest version of the RNAseq STAR-quantified data for all 787 cases performed on the GDC portal (https://gdc.cancer.gov/about-gdc/ contributed-genomic-data-cancer-research/foundation-medicine/multiple-myeloma-researchfoundation-mmrf, accessed on 20 March 2023, release date: 29 March 2022, version 32). STAR-aligned read groups were quantified using the two-pass method that generated the final alignments onto the GRCh38.p0 reference genome to calculate the gene-level RNAseq raw STAR-count data (i.e., STAR-aligned unstranded number of reads aligned per gene per sample), and fragments per kilobase of transcript per million mapped reads (FPKM) normalized to the upper quartile FPKM (FPKM-UQ). Gene expression profiles were downloaded from the archived MMRF CoMMpass dataset, as previously reported [25].
Patient-level clinical data were processed via functional tools provided in GenomicDat-aCommons_1. 18.0 and the case IDs were matched with RNAseq unique identifiers utilizing the metadata ("metadata.cart.2022-09-04.json") from the GDC portal converted into R data by running the utilities rjson_0.2.21 and stringr_1.4.0. The database included 766 ISS-staged patients, including 267 stage I, 276 stage II, and 223 stage III patients. We compared the mRNA expression levels for EGFR/ERBB1, ERBB2/HER2, and ERBB3 in each subset using the DESeq2 package (DESeq2_1.34.0 implemented using R version 4.1.2; R Foundation for Statistical Computing, Vienna, Austria. (1 November 2021)) [70]. Gene-level normalized, log 2 -transformed STAR-counts were calculated utilizing a generalized linear model that fits STAR-counts to negative binomial distribution to approximate the empirical distribution of the count data [25,70]. Statistical significance was assessed by testing the null hypothesis that there is no differential expression across the two sample groups (Log 2 fold change = 0) using the Wald test [25,70] reporting the test statistic and p-value for each gene. Genes were considered differentially expressed if the p-values were less than 0.05 after adjusting for multiple comparisons using the Benjamini and Hochberg method [25,70,71]. To visualize the gene expression levels in box plots, the normalized log 2 values were calculated from the RNAseq count data using the variance stabilization method supplied by the algorithms in vsn_3.62.0 [72]. Pairwise t-tests were performed to compare the expression of ERBB1 vs. ERBB2 and ERBB3 vs. ERBB2 across 766 patients and corrected for multiple comparisons using the Benjamini and Hochberg method.
We correlated the mRNA expression levels of ERBB2/HER2 in malignant plasma cells from 787 MM patients with the mRNA expression levels of 14 transcription factors known to activate ERBB2 expression [26][27][28][29][30][31][32][33][34]  Pairwise correlation coefficients were determined for 240 gene combinations using the FPKM-UQ data for each gene. Correlation coefficients were visualized on a heatmap color coded for positive correlations (red = +1) to negative correlations (blue = −1). The clustering algorithm identified co-regulated sets of mRNA expression values using the statistical package ggcorrplot_0.1.3 implemented in R software (R Foundation for Statistical Computing, Vienna, Austria). Significant correlations were identified for t-test p-values less than 0.05 and false discovery rates (FDRs) less than 0.05.

Analysis of Patient Outcomes According to ERBB2/HER2 mRNA Expression Levels
RNAseq and overall survival (OS) data were available for 787 MM patients. Progressionfree survival (PFS) information was available for 669 MM patients. The Kaplan-Meier (KM) method was employed utilizing the software packages survival_3.2-13, survminer_0.4.9, and survMisc_0.5.5 operated in the R environment to compare the progression-free survival (PFS) (n = 669), overall survival (OS) (n = 787), and cancer-related mortality (CRM) (n = 716) outcomes of patient subsets, as previously described [25]. Log-rank p-values less than 0.05 were deemed significant. Graphical representations of the survival curves were visualized using graph-drawing packages implemented in the R programming environment: dplyr_1.0.7, ggplot2_3.3.5, and ggthemes_4.2.4. The percentiles of patients expressing ERBB2/HER2 were determined using the metric fragments per kilobase of transcript per million mapped reads (FPKM) calculation normalized to the upper quartile FPKM (FPKM-UQ) for the whole geneset in each sample. PFS, OS, and cancer-related mortality (CRM) curves were compared for patients expressing high levels of ERBB2 (top 50th percentile of FPKM-UQ values) versus low levels of ERBB2 (bottom 50th percentile of FPKM-UQ values). Univariate Cox regression analyses were also performed to evaluate the impact of high-level ERBB2 expression on PFS and OS outcomes, as previously described [25]. We also investigated cancer-related mortality of four groups of patients with combinations of expression levels of ERBB1 and ERBB2: patients expressing high levels of both ERBB1 and ERBB2 (ERBB1 high /ERBB2 high ; top 50th percentile for both ERBB1 and ERBB2 mRNA expression); low levels of ERBB1 and ERBB2 (ERBB1 low /ERBB2 low ; below the 50th percentile of the lowest observed expression levels for both ERBB1 and ERBB2); and crossed combinations of ERBB1 and ERBB2 mRNA expression levels (ERBB1 low /ERBB2 high , and ERBB1 high /ERBB2 low ).
Multivariate analyses of the potential effect of the MMRF CoMMpass ERBB2/HER2 cohort groups on the PFS, OS, and cancer-related mortality outcomes were performed using the multivariate Cox proportional hazards model to adjust for other patient risk factors and staging criteria, as previously described [25]. Briefly, the model included (i) the mRNA expression level for ERBB2 as a linear co-variate or as a categorical variable comparing high versus low ERBB2 mRNA expression levels (50% cut-off for the range of FPKM-UQ values), (ii) the prognostic staging data according to the International Staging System (ISS), (iii) age, (iv) serum concentrations for albumin concentration and beta 2 microglobulin level, implemented in R (survival_3.2-13 ran in R version 4.1.2. Forest Plots were utilized, as previously reported [25]), to visualize the hazard ratios for Cox proportional hazards models for PFS, OS, and cancer-related mortality (survminer_0.4.9 ran in R version 4.1.2 (1 November 2021)). The life table hazard ratios (HRs) were estimated using the exponentiated regression coefficient for Cox proportional hazards analyses implemented in R (survival_3.2-13 ran in R version 4.1.2.).

Identifying Prognostically Relevant Signaling Proteins Networked to ERBB2 Expression Using the STRING Interaction Algorithm
Putative protein-protein interaction networks derived from mRNA levels with prognostic impact on the OS were constructed using the STRING11.5 algorithm (http://string-db.org/, accessed 20 May 2023) to identify candidate hub proteins connecting ERBB2 gene expression to potential networked interactions involving intracellular signaling molecules [71]. In these diagrams, the nodes depicted protein identifiers, and the edges depicted associations between the proteins.
Networks were visualized by first seeding the inputs and then growing the networks to identify connecting hubs between the inputs. The edges indicated the confidence level for the association calculated from various lines of experimental evidence (experiments, databases, and co-expression). Interaction scores greater than 0.7 were considered to define associations between the proteins and indicated by the thickness of the connecting lines for scores of 0.7 and 0.9. In the first instance, ERBB2 was seeded with the genes whose mRNA expression level was equal to or greater than the 50th percentile resulting in worse overall survival outcomes in 787 MM patients, namely, SHC1, CBL, PTPN11, FYN, KRAS, HRAS, NRAS, and MAP2K1.
These genes served as inputs to identify additional interactants by interrogating the STRING database for associations with scores equal to or greater than 0.7, no more than 20 interactants in the first shell, and no more than 5 interactants in the second shell of networked proteins. This iterative process identified a full module of 16 genes including ERBB2 whose augmented expression was highly associated with poor OS outcomes (p < 0.05, FDR < 0.1). A Kaplan-Meier overall survival analysis was performed to compare the top 50th percentile expression of each gene coding for the networked protein versus the bottom 50th percentile expression that resulted in 16 genes exhibiting worse survival outcomes (Log-rank p-values < 0.05, FDR = 0.06) at higher levels of expression in MM patients, namely: CBL, CRKL, EIF4EBP1, ERBB2, FGR, FKBP1A, FYN, HRAS, KRAS, MAP2K1, MTOR, NRAS, PTPN11, RALA, RPS6KB1, and SHC1. A Markov cluster (MCL) algorithm was applied to visualize subnetworks of highly connected nodes in the 16 genes coding for proteins shown to result in worse survival outcomes in MM patients. The granularity of visualizing the subnetworks is controlled via the Markov cluster inflation factor ranging from 1.1 to 10 where, in our analysis, we discerned 2 highly connected sub-networks by setting the value to 2.1 for the 16 genes examined.

Conclusions
In MM patients, ERBB2 was expressed at significantly higher levels than ERBB1 as well as ERBB3 across all three stages of the disease. Upregulated expression of ERBB2 mRNA in MM cells was correlated with amplified expression of mRNAs for transcription factors that recognize the ERBB2 gene promoter sites. Patients with higher levels of ERBB2 mRNA in their malignant plasma cells experienced significantly increased cancer mortality, shorter progression-free survival, and worse overall survival than other patients. The adverse impact of high ERBB2 expression on patient survival outcomes remained significant in multivariate Cox proportional hazards models that accounted for the effects of other prognostic factors. Funding: This study received funding from Ares Pharmaceuticals. Authors F.M.U. and S.Q. participated in the analysis and decision to submit the manuscript for publication are affiliated with the funder Ares Pharmaceuticals.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All data are presented in the manuscript. The publicly available archived databases that were used to generate the data are provided in Section 2. Gene level RNAseq raw count data (STAR-aligned unstranded number of reads aligned per gene per sample) were downloaded from the archived MMRF CoMMpass dataset by connecting to the GDC portal (https://gdc.cancer.gov/about-gdc/contributed-genomic-data-cancer-research/foundationmedicine/multiple-myeloma-research-foundation-mmrf (accessed on 24 September 2022)) using the Bioconductor package, GenomicDataCommons_1.18.0 (release date 11 August 2021 with full functionality as provided by TCGAbiolinks for accessing GDC data (https://bioconductor.org/packages/ release/bioc/html/GenomicDataCommons.html (accessed on 20 September 2022)) implemented in R version 4.1.2 (1 November 2021). The mRNA expression data were deposited in files appended with " . . . .rna_seq.augmented_star_gene_counts.tsv". Clinical data for each MM patient were also acquired via functions provided in GenomicDataCommons_1. 18.0 and the case IDs were matched with RNAseq unique identifiers utilizing the metadata ("metadata.cart.2022-09-04.json") from the GDC portal converted into R data by running the utilities rjson_0.2.21 and stringr_1.4.0.

Conflicts of Interest:
The author F.M.U. is employed by Ares Pharmaceuticals, and the author S.Q. serves as a consultant for Ares Pharmaceuticals. All authors declare no other competing interest.